선형 회귀 분석의 기초

  • 결정론적 모형은 그냥 함수를 찾는 것. 간단한 함수부터 시작을 한다. 간단한 함수는 선형식을 의미하는 듯
  • 선형 회귀 분석은 부호, 크기, 관계 등을 알려주기 때문에 불안전하다는 단점에도 불구하고 잘 쓰이고 있다. 비선형회귀분석의 문제점으로는 overfitting 현상이 발생한다는 점. 그리고 방법도 너무 많다는 점
  • cross validation이란 x를 남겨두는 것. 함수 검사를 위해서. 진짜 시험을 남겨두는 것과 같은 원리. 적어도 3개 이상 남겨둔다.

회귀 분석(regression analysis)은 입력 자료(독립 변수) $x$와 이에 대응하는 출력 자료(종속 변수) $y$간의 관계를 정량화 하기 위한 작업이다.

회귀 분석에는 결정론적 모형(Deterministic Model)과 확률적 모형(Probabilistic Model)이 있다.

결정론적 모형은 단순히 독립 변수 $x$에 대해 대응하는 종속 변수 $y$를 계산하는 함수를 만드는 과정이다.

$$ \hat{y} = f \left( x; \{ x_1, y_1, x_2, y_2, \cdots, x_N, y_N \} \right) = f (x; D) = f(x) $$

여기에서 $ \{ x_1, y_1, x_2, y_2, \cdots, x_N, y_N \} $ 는 모형 계수 추정을 위한 과거 자료이다.

만약 함수가 선형 함수이면 선형 회귀 분석(linear regression analysis)이라고 한다.

$$ \hat{y} = w_0 + w_1 x_1 + w_2 x_2 + \cdots + w_D x_D $$

Augmentation(증가 개념)

일반적으로 회귀 분석에 앞서 다음과 같이 상수항을 독립 변수에 포함하는 작업이 필요할 수 있다. 이를 feature augmentation이라고 한다.

$$ x_i = \begin{bmatrix} x_{i1} \\ x_{i2} \\ \vdots \\ x_{iD} \end{bmatrix} \rightarrow x_{i,a} = \begin{bmatrix} 1 \\ x_{i1} \\ x_{i2} \\ \vdots \\ x_{iD} \end{bmatrix} $$

augmentation을 하게 되면 모든 원소가 1인 벡터를 feature matrix 에 추가된다.

$$ X = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1D} \\ x_{21} & x_{22} & \cdots & x_{2D} \\ \vdots & \vdots & \vdots & \vdots \\ x_{N1} & x_{N2} & \cdots & x_{ND} \\ \end{bmatrix} \rightarrow X_a = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1D} \\ 1 & x_{21} & x_{22} & \cdots & x_{2D} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_{N1} & x_{N2} & \cdots & x_{ND} \\ \end{bmatrix} $$

augmentation을 하면 가중치 벡터(weight vector)도 차원이 증가하여 전체 수식이 다음과 같이 단순화 된다.

$$ w_0 + w_1 x_1 + w_2 x_2 = \begin{bmatrix} 1 & x_1 & x_2 \end{bmatrix} \begin{bmatrix} w_0 \\ w_1 \\ w_2 \end{bmatrix} = x_a^T w $$

In [1]:
from sklearn.datasets import make_regression
bias = 100
X0, y, coef = make_regression(n_samples=100, n_features=1, bias=bias, noise=10, coef=True, random_state=1)
X = np.hstack([np.ones_like(X0), X0])
X[:5]


Out[1]:
array([[ 1.        , -0.61175641],
       [ 1.        , -0.24937038],
       [ 1.        ,  0.48851815],
       [ 1.        ,  0.76201118],
       [ 1.        ,  1.51981682]])

OLS (Ordinary Least Squares)

OLS는 가장 기본적인 결정론적 회귀 방법으로 Residual Sum of Squares(RSS)를 최소화하는 가중치 벡터 값을 미분을 통해 구한다.

  • Residual 잔차 $$ e_i = {y}_i - x_i^T w $$
  • Stacking (Vector Form) $$ e = {y} - Xw $$
  • Residual Sum of Squares (RSS) $$\begin{eqnarray} \text{RSS} &=& \sum (y_i - \hat{y}_i)^2 \\ &=& \sum e_i^2 = e^Te \\ &=& (y - Xw)^T(y - Xw) \\ &=& y^Ty - 2y^T X w + w^TX^TXw \end{eqnarray}$$
  • Minimize using Gradient $$ \dfrac{\partial \text{RSS}}{\partial w} = -2 X^T y + 2 X^TX w = 0 $$
$$ X^TX w = X^T y $$$$ w = (X^TX)^{-1} X^T y $$

여기에서 그레디언트를 나타내는 다음 식을 Normal equation 이라고 한다.

$$ X^T y - X^TX w = 0 $$

Normal equation 에서 잔차에 대한 다음 특성을 알 수 있다.

$$ X^T (y - X w ) = X^T e = 0 $$
  • bias는 상수항이자 y절편. 사이킷은 내부에 있고 stats모델에서는 명령어 하나 불러야 한다?

In [2]:
y = y.reshape(len(y), 1)
w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y)
print("bias:", bias)
print("coef:", coef)
print("w:\n", w)


bias: 100
coef: 80.71051956187792
w:
 [[ 102.02701439]
 [  81.59750943]]

In [4]:
w = np.linalg.lstsq(X, y)[0]
w


Out[4]:
array([[ 102.02701439],
       [  81.59750943]])

In [5]:
xx = np.linspace(np.min(X0) - 1, np.max(X0) + 1, 1000)
XX = np.vstack([np.ones(xx.shape[0]), xx.T]).T
yy = np.dot(XX, w)
plt.scatter(X0, y)
plt.plot(xx, yy, 'r-')
plt.show()


scikit-learn 패키지를 사용한 선형 회귀 분석

sklearn 패키지를 사용하여 선형 회귀 분석을 하는 경우에는 linear_model 서브 패키지의 LinearRegression 클래스를 사용한다.

  • 입력 인수

    • fit_intercept : 불리언, 옵션
      • 상수상 추가 여부
    • normalize : 불리언, 옵션
      • 회귀 분석전에 정규화 여부
  • 속성

    • coef_ : 추정된 가중치 벡터
    • intercept_ : 추정된 상수항

Diabetes Regression


In [1]:
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
dfX_diabetes = pd.DataFrame(diabetes.data, columns=["X%d" % (i+1) for i in range(np.shape(diabetes.data)[1])])
dfy_diabetes = pd.DataFrame(diabetes.target, columns=["target"])
df_diabetes0 = pd.concat([dfX_diabetes, dfy_diabetes], axis=1)
df_diabetes0.tail(3)


Out[1]:
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 target
439 0.041708 0.050680 -0.015906 0.017282 -0.037344 -0.013840 -0.024993 -0.011080 -0.046879 0.015491 132.0
440 -0.045472 -0.044642 0.039062 0.001215 0.016318 0.015283 -0.028674 0.026560 0.044528 -0.025930 220.0
441 -0.045472 -0.044642 -0.073030 -0.081414 0.083740 0.027809 0.173816 -0.039493 -0.004220 0.003064 57.0

In [2]:
from sklearn.linear_model import LinearRegression
model_diabets = LinearRegression().fit(diabetes.data, diabetes.target)
print(model_diabets.coef_)
print(model_diabets.intercept_)


[ -10.01219782 -239.81908937  519.83978679  324.39042769 -792.18416163
  476.74583782  101.04457032  177.06417623  751.27932109   67.62538639]
152.133484163

In [3]:
predictions = model_diabets.predict(diabetes.data)
plt.scatter(diabetes.target, predictions)
plt.xlabel("prediction")
plt.ylabel("target")
plt.show()



In [4]:
mean_abs_error = (np.abs(((diabetes.target - predictions) / diabetes.target)*100)).mean()
print("MAE: %.2f%%" % (mean_abs_error))


MAE: 38.79%

In [5]:
sk.metrics.median_absolute_error(diabetes.target, predictions)


Out[5]:
38.522833662660567

In [6]:
sk.metrics.mean_squared_error(diabetes.target, predictions)


Out[6]:
2859.6903987680657

Boston Housing Price


In [7]:
from sklearn.datasets import load_boston
boston = load_boston()
dfX_boston = pd.DataFrame(boston.data, columns=boston.feature_names)
dfy_boston = pd.DataFrame(boston.target, columns=["MEDV"])
df_boston0 = pd.concat([dfX_boston, dfy_boston], axis=1)
df_boston0.tail(3)


Out[7]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT MEDV
503 0.06076 0.0 11.93 0.0 0.573 6.976 91.0 2.1675 1.0 273.0 21.0 396.90 5.64 23.9
504 0.10959 0.0 11.93 0.0 0.573 6.794 89.3 2.3889 1.0 273.0 21.0 393.45 6.48 22.0
505 0.04741 0.0 11.93 0.0 0.573 6.030 80.8 2.5050 1.0 273.0 21.0 396.90 7.88 11.9

In [8]:
model_boston = LinearRegression().fit(boston.data, boston.target)
print(model_boston.coef_)
print(model_boston.intercept_)


[ -1.07170557e-01   4.63952195e-02   2.08602395e-02   2.68856140e+00
  -1.77957587e+01   3.80475246e+00   7.51061703e-04  -1.47575880e+00
   3.05655038e-01  -1.23293463e-02  -9.53463555e-01   9.39251272e-03
  -5.25466633e-01]
36.4911032804

In [9]:
predictions = model_boston.predict(boston.data)
plt.scatter(predictions, boston.target)
plt.xlabel("prediction")
plt.ylabel("target")
plt.show()



In [10]:
mean_abs_error = (np.abs(((boston.target - predictions) / boston.target)*100)).mean()
print("MAE: %.2f%%" % (mean_abs_error))


MAE: 16.43%

In [11]:
sk.metrics.median_absolute_error(boston.target, predictions)


Out[11]:
2.455946432803632

In [12]:
sk.metrics.mean_squared_error(boston.target, predictions)


Out[12]:
21.897779217687503

statsmodels 를 사용한 선형 회귀 분석

  • 실제로는 선형회귀분석의 경우 이 모델을 사용한다.

statsmodels 패키지에서는 OLS 클래스를 사용하여 선형 회귀 분석을 실시한다.

statsmodels.regression.linear_model.OLS(endog, exog=None)

  • 입력 인수
    • endog : 종속 변수. 1차원 배열
    • exog : 독립 변수, 2차원 배열.

statsmodels 의 OLS 클래스는 자동으로 상수항을 만들어주지 않기 때문에 사용자가 add_constant 명령으로 상수항을 추가해야 한다.

모형 객체가 생성되면 fit, predict 메서드를 사용하여 추정 및 예측을 실시한다.

예측 결과는 RegressionResults 클래스 객체로 출력되면 summary 메서드로 결과 보고서를 볼 수 있다.


In [14]:
df_diabetes = sm.add_constant(df_diabetes0)
df_diabetes.tail(3)


Out[14]:
const X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 target
439 1 0.041708 0.050680 -0.015906 0.017282 -0.037344 -0.013840 -0.024993 -0.011080 -0.046879 0.015491 132.0
440 1 -0.045472 -0.044642 0.039062 0.001215 0.016318 0.015283 -0.028674 0.026560 0.044528 -0.025930 220.0
441 1 -0.045472 -0.044642 -0.073030 -0.081414 0.083740 0.027809 0.173816 -0.039493 -0.004220 0.003064 57.0

In [15]:
model_diabets2 = sm.OLS(df_diabetes.ix[:, -1], df_diabetes.ix[:, :-1])
result_diabetes2 = model_diabets2.fit()
result_diabetes2


Out[15]:
<statsmodels.regression.linear_model.RegressionResultsWrapper at 0xb3bf160>
  • DF는 ndarray의 리스트
    • 리스트 of 벡터
    • ndarray는 다 된다
    • list of list는 리스트 안에 리스트
  • Dep. Variable은 우리가 구할 값
    • target은 라벨
  • No. Observations은 샘플수
  • Df Model은 parameter-1
  • std err는 coef의 +- err 수치
  • 가장 먼저 P>|t| 이거부터. 이게 중요. 0인지 아닌지. 0이면 살리고 아니면 죽일 가능성이 높다.
  • Prob(Omnibus) = 0.471이면 그냥 정규분포다.
  • Cond. No가 10000 이하면 괜찮아
  • 그 다음 봐야 할 것은 coef. coef_ : 추정된 가중치 벡터. ‘-’이면 악영향
  • Result가 별도로 저장되는 것은 stats 모델의 특징. 사이킷은 아니야

In [16]:
print(result_diabetes2.summary())


                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 target   R-squared:                       0.518
Model:                            OLS   Adj. R-squared:                  0.507
Method:                 Least Squares   F-statistic:                     46.27
Date:                Wed, 07 Sep 2016   Prob (F-statistic):           3.83e-62
Time:                        08:12:22   Log-Likelihood:                -2386.0
No. Observations:                 442   AIC:                             4794.
Df Residuals:                     431   BIC:                             4839.
Df Model:                          10                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const        152.1335      2.576     59.061      0.000       147.071   157.196
X1           -10.0122     59.749     -0.168      0.867      -127.448   107.424
X2          -239.8191     61.222     -3.917      0.000      -360.151  -119.488
X3           519.8398     66.534      7.813      0.000       389.069   650.610
X4           324.3904     65.422      4.958      0.000       195.805   452.976
X5          -792.1842    416.684     -1.901      0.058     -1611.169    26.801
X6           476.7458    339.035      1.406      0.160      -189.621  1143.113
X7           101.0446    212.533      0.475      0.635      -316.685   518.774
X8           177.0642    161.476      1.097      0.273      -140.313   494.442
X9           751.2793    171.902      4.370      0.000       413.409  1089.150
X10           67.6254     65.984      1.025      0.306       -62.065   197.316
==============================================================================
Omnibus:                        1.506   Durbin-Watson:                   2.029
Prob(Omnibus):                  0.471   Jarque-Bera (JB):                1.404
Skew:                           0.017   Prob(JB):                        0.496
Kurtosis:                       2.726   Cond. No.                         227.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

In [17]:
df_boston = sm.add_constant(df_boston0)
model_boston2 = sm.OLS(df_boston.ix[:, -1], df_boston.ix[:, :-1])
result_boston2 = model_boston2.fit()
print(result_boston2.summary())


                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   MEDV   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Wed, 07 Sep 2016   Prob (F-statistic):          6.95e-135
Time:                        08:20:04   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         36.4911      5.104      7.149      0.000        26.462    46.520
CRIM          -0.1072      0.033     -3.276      0.001        -0.171    -0.043
ZN             0.0464      0.014      3.380      0.001         0.019     0.073
INDUS          0.0209      0.061      0.339      0.735        -0.100     0.142
CHAS           2.6886      0.862      3.120      0.002         0.996     4.381
NOX          -17.7958      3.821     -4.658      0.000       -25.302   -10.289
RM             3.8048      0.418      9.102      0.000         2.983     4.626
AGE            0.0008      0.013      0.057      0.955        -0.025     0.027
DIS           -1.4758      0.199     -7.398      0.000        -1.868    -1.084
RAD            0.3057      0.066      4.608      0.000         0.175     0.436
TAX           -0.0123      0.004     -3.278      0.001        -0.020    -0.005
PTRATIO       -0.9535      0.131     -7.287      0.000        -1.211    -0.696
B              0.0094      0.003      3.500      0.001         0.004     0.015
LSTAT         -0.5255      0.051    -10.366      0.000        -0.625    -0.426
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                     1.51e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.51e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

RegressionResults 클래스는 분석 결과를 다양한 속성에 저장해주므로 추후 사용자가 선택하여 활용할 수 있다.


In [18]:
dir(result_boston2)


Out[18]:
['HC0_se',
 'HC1_se',
 'HC2_se',
 'HC3_se',
 '_HCCM',
 '__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_cache',
 '_data_attr',
 '_get_robustcov_results',
 '_is_nested',
 '_wexog_singular_values',
 'aic',
 'bic',
 'bse',
 'centered_tss',
 'compare_f_test',
 'compare_lm_test',
 'compare_lr_test',
 'condition_number',
 'conf_int',
 'conf_int_el',
 'cov_HC0',
 'cov_HC1',
 'cov_HC2',
 'cov_HC3',
 'cov_kwds',
 'cov_params',
 'cov_type',
 'df_model',
 'df_resid',
 'diagn',
 'eigenvals',
 'el_test',
 'ess',
 'f_pvalue',
 'f_test',
 'fittedvalues',
 'fvalue',
 'get_influence',
 'get_robustcov_results',
 'initialize',
 'k_constant',
 'llf',
 'load',
 'model',
 'mse_model',
 'mse_resid',
 'mse_total',
 'nobs',
 'normalized_cov_params',
 'outlier_test',
 'params',
 'predict',
 'pvalues',
 'remove_data',
 'resid',
 'resid_pearson',
 'rsquared',
 'rsquared_adj',
 'save',
 'scale',
 'ssr',
 'summary',
 'summary2',
 't_test',
 'tvalues',
 'uncentered_tss',
 'use_t',
 'wald_test',
 'wresid']

statsmodel는 다양한 회귀 분석 결과 플롯도 제공한다.

  • plot_fit(results, exog_idx) Plot fit against one regressor.
  • abline_plot([intercept, ...]) Plots a line given an intercept and slope.
  • influence_plot(results[, ...]) Plot of influence in regression.
  • plot_leverage_resid2(results) Plots leverage statistics vs.
  • plot_partregress(endog, ...) Plot partial regression for a single regressor.
  • plot_ccpr(results, exog_idx) Plot CCPR against one regressor.
  • plot_regress_exog(results, ...) Plot regression results against one regressor.

In [19]:
sm.graphics.plot_fit(result_boston2, "CRIM")
plt.show()